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Abstract 

The paper describes the different methods, used in the MAGIC experiment, to 
unfold experimental energy distributions of cosmic ray particles (7-rays). Ques- 
tions and problems related to the unfolding are discussed. Various procedures are 
proposed which can help to make the unfolding robust and reliable. The different 
methods and procedures are implemented in the MAGIC software and are used in 
most of the analyses. 
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1 Introduction 

In an Imaging-Air-Cherenkov- Telescope (lACT) experiment like MAGIC [1] 
the energy E of the cosmic ray particle (7-ray) is not exactly known. It has 
to be estimated, the energy resolution being in the order of 20 to 40%. As 
a consequence, the experimentally measured energy spectrum is biased. The 
procedure to correct for the effects due to the finite energy resolution is called 
unfolding. 

While in high-energy-physics experiments unfolding is a widely used technique, 
this is not the case in present day's lACT experiments. This paper deals with 
the unfolding procedure, which is applied as a standard tool in the MAGIC 
experiment. The different unfolding methods are explained in detail. Emphasis 
is put on the discussion of questions and problems related to the application 
of the unfolding to real data. It is not the aim of the paper to give a complete 
derivation of all formulas. For this the reader is referred to the publications [2] 
to [18]. An excellent review of unfolding methods is given in [13]. The present 
paper makes use of many ideas discussed in that paper. 

Although here only differential energy spectra are considered, the procedures 
are equally well applicable to distributions of other quantities, including dis- 
tributions in more than one dimension [19]. 

The layout of the paper is as follows. In Section 2 the notation is defined and 
the motivation for the unfolding procedure is given. The different unfolding 
methods, which means the different ways of regularization, are presented in 
Section 3. The so-called Forward Unfolding, which represents an imphcit un- 
folding under the assumption of a certain parametrization of the unfolding 
result is explained in Section 4. In Section 5 two quantities are introduced, 
which are useful for an optimal choice of the regularization strength. The 
criteria for this choice are collected in Section 6. Section 7 discusses various 
technical aspects which are important in the application of the unfolding pro- 
cedure to real data. Two particular technical procedures, which ensure an 
unbiased and robust unfolding, are presented in Sections 8 and 9. Some un- 
folding results, obtained by applying the unfolding procedure to data taken 
in the MAGIC experiment, are discussed in Section 10. Finally the pros and 
contras of the method of Correction Factors, which is an alternative way of 
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correcting data for effects due to the finite experimental resolution, are listed 
in Section 11. A summary is given in Section 12. 



2 The aim of the unfolding procedure 

In this Section the notations are defined and the motivation for the unfolding 
procedure is explained. 



2. 1 Notation 

The true and measured (estimated) values of the energy of the cosmic ray 
particle are denoted by Et^ue and E^st respectively. The data are assumed to 
be binned in histograms, and certain binnings are chosen independently for 
the distributions in Etme and E^st- Furthermore, the following definitions are 
introduced: 





number of events in bin i of Ee,st 
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with 

na 

^M,,- = 1 for all J (1) 

i=l 

In the following it is assumed that the rank nr of G is equal to (and not less 
than) the minimum of na and nb, where na and nb are the number of bins in 
Egst and Efrue respectively, which are used in the unfolding. This can always 
be achieved by a proper choice of the binnings in Etme and Eest- 

The migration matrix contains the most likely fraction of events moving from 
a bin j in Efme into a bin i of Eggf, due to the finite experimental energy 
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resolution: 

nb 

Yi = Y^Mij-Sj (i = l,...na) 
or in matrix notation Y — M ■ S (2) 

The migration matrix M is obtained from Monte Carlo (MC) simulations, 
in which the development of the air shower (induced by the cosmic ray par- 
ticle in the atmosphere), the emission of Cherenkov hght in the air shower, 
the geometrical, optical and electronic properties of the telescope and the ex- 
perimental procedures in the data analysis (shower reconstruction, 7/hadron 
separation, energy estimation, selections, cuts) are simulated [20]. M is com- 
puted from a 2-dimensional plot of the number of reconstructed MC events 
in the Eest-Etme plane, which was produced under the same conditions (selec- 
tions, cuts) as the experimental distribution Y . It is the aim of the unfolding 
procedure to determine the true distribution given Y and M. 

It should be noted that M only describes the migration of events. It does not 
describe losses of events, which will occur due to the finite acceptance of the 
detector, due to the trigger conditions and due to additional selections and 
cuts. In an lACT experiment these losses are also determined by Monte Carlo 
simulations, and the corresponding correction factor is the effective collection 
area A^ff^Etme)- ^e//(-£'true) has to be computed again under the same condi- 
tions as the experimental distribution Y . Apart from minor effects (see Section 
8) , this correction can be performed quite independently of the unfolding. 

The unfolding can be understood as a reshuffling of events from the bins of E^st 
into the bins of Etme- In this procedure the numerical values of the bin edges, 
both for Etrue and £Jest, are completely irrelevant. Etme and Eest may be even 
two different physical quantities, with completely different ranges of values 
and different units, like Et^ue = (the true energy of the cosmic ray particle) 
and Egst = (the total number of Cherenkov photons measured for the shower) 
[21]. Of course, unfolding makes only sense if Etrue and E^st are sufficiently 
strongly correlated, otherwise the distribution of Etme cannot be inferred from 
a distribution of Eest- This is in contrast to the method of Correction Factors 
(see Section 11), where E^gt has to be a good estimate of Etme in any case. 

Because the binnings in Etme and Eggt can be chosen independently and to a 

certain degree arbitrarily (see Section 7), the migration matrix M is in general 
not square. For this reason eq. (2) can in general not be inverted to obtain S 
as • Y. 
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As the unfolding is equivalent to a reshuffling of events from the bins of the 
measured distribution into the bins of the true distribution the unfolding is 
not restricted to 1-dimensional distributions but can in the same way be ap- 
plied to multi-dimensional distributions. The information necessary for the 
unfolding procedure is completely contained in the corresponding migration 
matrix, which in the case of multi-dimensional distributions describes the mi- 
gration of events from the bins of the true multi-dimensional distribution into 
the bins of the measured multi- dimensional distribution. The dimensions of 
the measured and true distributions may also be different. An example for an 
unfolding in 2 dimensions is given in [19]. 



2.2 The direct solution of Y = M ■ S 



Very generally, the solution S of the system of linear equations (2) can be 
obtained by minimizing the Least-Squares expression 



Xo = {Y-M-Sf-K-^-{Y-M-S) (3) 

where the nb components of S are the free parameters. Minimizing Xo will 
yield solutions for S which, after folding with M, are best compatible with 
the measurement Y. 



Two cases have to be distinguished: 

• The underconstrained case nr = na < nb. 

Because nr = na the na x na matrix G can be inverted and a particular 
solution 5*0 can be written as 

^M^ -C with C ^ G-^-Y (4) 

If na < nb, the solutions S = Sq + St form a space of {nb — na) dimensions 
with M ■ St — 0. For na — nb eq.(4) reduces to 

5° = -G-^ -Y = M-^ ■ Y (5) 

and S^ is the only solution. In both cases M-S — M-S'^ — Y, implying xl ~ 
0. The solutions are independent of the covariance matrix K. Moreover, 
because of (1) the total number of events is not changed: 

= EE^^.-^? = E^^ (6) 

i i j 3 

• The overconstrained case na > nb = nr. 
Minimizing Xq by varying S yields 
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(7) 



where H is the nb x nb matrix 



H = M^K-^M 



(8) 



The solution S^^^ now depends on K. The minimum value of Xq becomes 



Expression (7) is also valid for nr = na = nb, in which case it reduces to 
gLSQ ^ M-^-Y. 

It can be shown that the direct solutions (4) and (7) may lead to large errors 
of Sj, reflected in large absolute values of the elements of the error matrix T 
of S. This behaviour can be traced back to small eigenvalues of the matrix G 
and H respectively [13]. 



3 Unfolding with Regular izat ion - The different Unfolding Meth- 



In order to reduce the large errors of S, a procedure called regularization is ap- 
plied. By the regularization additional constraints are imposed on S, by which 
some information in the measurements Yi is discarded. Regularization can be 
viewed as a smearing of the unfolded distribution with some finite resolution, 
which reduces the correlations between the Si of adjacent bins at the expense 
that Si is no longer an unbiased estimate of the true distribution [17]. 

The bias increases with increasing regularization strength. Nevertheless, it 
turns out that with properly tuned regularization (sec Section 6) solutions 
can be obtained which are much closer to the true distribution than the direct 
solutions (4) or (7). 

It is evident that regularization is particularly important in the undercon- 
strained case. However, also in the overconstrained case regularization makes 
sense: even if the system of equations (2) is formally overconstrained {na > 
nb), it may be effectively underconstrained. This happens for example if some 
of the measurements Yi have much larger errors than the other Yi. 

In the following, three different ways of regularization are described [13] . 



Xl ^ {Y-M- S'^^'^f -K-^-iY-M- 



(9) 



ods 
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3. 1 Adding a regularization term in the expression for Xo 



In some unfolding methods regularization is performed by adding a regular- 
ization term Reg{S) in the expression for Xq (eq.3) 



= 2 - x^ + Regis) (10) 

also called regularization parameter, is a weight which allows to steer the 
regularization strength: large values of w correspond to weak regularization, 
small values to strong regularization. 



• Tikhonov's method 

In Tikhonov's method [4] the regularization term is defined as 

R^9{s) - Y.[-^] (11) 



(d^S\ 

For the second derivative I ^ ^ I of 5 in bin j different approximations 
may be used. The expression used in the MAGIC software [22] is 

This is actually an approximation for the bin-to- bin variation of AS/S. 

( d '^S\ 

In [7] ( -^-Y j is calculated from a sphne representation of S. 

For a given value of w and after specifying Reg{S), expression (10) can be 
minimized numerically by varying the components of S. The minimization 
also provides the error matrix T of S. The regularization matrix R, defined 
in (24), can be calculated numerically by performing minimizations with 
modified values of Yi. 

Schmelling's method 

In this method, which is discussed in great detail in [14, 17], the regular- 
ization term is set equal to the "cross entropy" 

nb 

Regis) = EP.-ln^ (13) 

j=i 

Pj is the normalized distribution S 
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and e is a normalized prior distribution, which describes a prior knowledge 
about S. The cross entropy Reg{S) quantifies by how much p deviates from 
e. Finding S by minimizing the cross entropy Reg{S) simultaneously with 
the least squares expression Xo called the method of " Reduced Cross En- 
tropy" . 

With Reg{S) from (13), expression (10) can now be minimized to ob- 
tain the unfolded distribution S. Note that all components of S and of the 
prior distribution are required to be > 0, because otherwise Reg{S) in (13) 
cannot be defined. The expressions for the error matrix T oi S and for the 
regularization matrix R are given in [14, 17]. 

In the MAGIC software, the condition J2i = J2j Sj is used as an additional 
constraint, when minimizing x^- By this one degree of freedom is gained. 

3.2 Spectral Window method 

In some unfolding methods regularization is performed by suppresing small 
eigenvalues A/ of G by a factor /(A/) [11]. By the suppression factor f{Xi) the 
matrix G and its inverse are modified. In terms of the eigenvectors gi of G 
they read 

nr 

G = Y. fih) ■ Xi ■ gigj (15) 

I 

G^l^^m.g^gT (16) 

I A; 

where the sums extend ovcivall eigenvalues A; which are different from zero. 
Like G and G^^, G and G~^ are na x na matrices. Without suppression, 
/(A;) = 1, G~^ is equal to G~^ in the underconstrained case nr = na < nb. 
In the overconstrained case, na > nb — nr, G~^ is undefined but G~^ can be 
calculated. 

A similar factor /(/t/) can be defined to suppress small eigenvalues ki of H (eq. 
8). There is considerable freedom as to the choice of the values or expressions 
for /(A;) and /(«;). One may introduce a parameter i such that in the limit 
i — >■ oo the suppression factors f{Xi,i) and f{Ki,i) tend to 1. i has a similar 
meaning as the weight w in eq.(lO): it determines the regularization strength 
and for i — oo the solutions tend to the direct solutions (4) and (7) respec- 
tively. 
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The expressions for the error matrix T oi S and for the regularization matrix 
R are given in [13]. 



3.3 Regularization by iteration 



Another way of unfolding is to calculate a solution S iteratively [3, 9, 11, 12]. 
The regularization is done by stopping the iteration at some point. In this 
case the number of iterations i plays a similar role as the weight w in (10). 
In the limit of an infinite number of iterations, which is equivalent to a very 
large weight w, the solution tends to the direct solution (4) or (7) respectively. 



nr = na < nb 

A simple iteration scheme [13] for solving Y — M ■ S with -Y — C is 



C 



C' - T- {GC' - Y) 



(17) 



where i is the iteration number and r is a relaxation parameter. The latter 
should be chosen in the range < r < 2/\max-i where Xmax is the largest 
eigenvalue of G. Eq.(17) leads to 



i-l 



C' = (1 - TGf • C° + r • ^(1 - TGy ■ Y 

3=0 



(18) 



where is the starting value of C. The unfolded distribution is obtained 

as 5' = M^-C\ 



In terms of the suppression factor f{Xi,i) this regularization can be ex- 
pressed as [13] 



/(Am) 



i-a-rx. 



if C° is set to zero, and 

f{Xi,i) = [l -(1-tAO' + {l-rXiYXi 



(19) 



(20) 



if C° is chosen to be equal to Y. The solutions 5** tends to the direct solution 
S^'^M'^ -C^M"^ -G-^-Y (eq. (4)). 

This procedure can also be applied in the overconstrained case, na > 
nb = nr. In this case G~^ is undefined and^ has to be replaced by G~^ 
(eq.l6). The solutions S' tend to S = ■ G-^ ■ F, with /(A;) = 1. 
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• na > nb = nr 

In a similar way one may define suppression factors /{ni, i) for tlie eigenval- 
ues of the matrix H. In this case the solution S tends to the solution S^^^ 
(eq. (7)). 

In the last 2 methods, Spectral Window method and Regularization by itera- 
tion, Xq doesn't appear explicitly. However, its value can be calculated and it 
is taken into account when determining the optimum regularization strength 
(see Section 6). 

Once the the unfolded distribution Sk is determined, by any of the methods 
described in this Section, the differential energy spectrum $fc of 7-rays is 
calculated using eq. (22). $fc has the meaning of the average differential 7-ray 
flux in the k-ih. bin of Etme (eq.(23)). 



4 Forward Unfolding 



An implicit unfolding can be done by representing 5 as a parametric function 
^kio.) = f{Ek', q) with parameters q — {qi, q2, ■■■qnq) and minimizing Xq in (3) 



na / nb \ / nb \ 

Xo = Y.[Y^-Y.M^k■ Su{q) ■ {k-') • - E Mji ■ Si{q) (21) 

i,j=l \ k=l / ■' \ 1=1 J 

with respect to the parameters q. The number of measurements is equal to na, 
the number of unknowns nq. Thus, the problem is overconstrained if nq < na, 
independent of the value of nb. One degree of freedom is gained if the total 
number of events is required to stay constant: J2j Sj — J2i'^i- 

In many cases the minimization of Xo can be performed analytically, by solving 

— — = 0, similarly to the procedure described in Section 2.2. 
oq 

The parametrization of S can be written in the form 



S,{q) = ■ Teff ■ A(ELe) • ^e// " ^Ldct " ^W^on (22) 

with 
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Here ^{Etme, q) is the assumed parametrization of the differential energy spec- 
trum of 7-rays, T^ff is the effective observation time, A{E^,^.^^,,) denotes the A;-th 
bin in Etrue or its width, ^g// is the effectice collection area, and A'^^^^^^ is the 
reduction factor due to an additional cut (for example in E^su see Section 7.3). 
A further correction A'l^g^p^^g^ can be introduced, if ^{Etme-i q) is supposed to 
represent the differential energy spectrum of 7-rays before absorption, either 
at the 7-source or by interaction with the extragalactic photon background. 
^addcut t)e determined from MC simulations, whereas A^f^g^^p^^^^^ can be 
calculated in models about the extragalactic photon background [23]. 

Parametrizing S as an analytic function of Efrue, with some free parameters 
q, can be understood as a kind of regularization, because it forces the solution 
S and its derivatives to be continuous, leading to a suppression of the noise 
component of S. 

The Forward Unfolding does not provide an unfolded distribution S. It pro- 
vides those parameter values q for an assumed parametrization of ^{Efrue), 
which minimize xl iii (21)- Of course, S can then be calculated from ^{Etrue) 
via (22). 

Under the assumption of a certain parametrization of ^{Efrue), the Forward 
Unfolding is a very robust method of determining the best parameter values q. 
Moreover, since there is no regularization strength to be adjusted, the uncer- 
tainty as to its choice does not exist. Therefore Forward Unfolding represents 
a powerful and useful check of the unfolding results obtained by any of the 
methods described in Section 3. In those methods the parametrization is only 
introduced after the actual unfolding of the measurements Y. 



5 Useful quantities in the unfolding 

In this Section two quantities are explained which are useful for judging the 
quality of the unfolding result: The error matrix (covariance matrix) T of the 
unfolded distribution S and the regularization matrix R. 

5. 1 The covariance matrix of the unfolded distribution S 

In all unfolding methods the covariance matrix T oi S can be determined. 
The trace of T, Trace{T), measures the noise component of -S", as Trace{K) 
measures the noise component of Y. 

In the methods where the solution S is given as S = D-Y, like in Schmelling's 
method or in the Spectral Window method, the covariance matrix T of 5" is 
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obtained by T = D-K-D^. 



In those methods where the solution is determined by a numerical minimiza- 
tion of x^, like in Tikhonov's method, T is obtained from the shape of in 
the region around the minimum. 



5.2 The regularization matrix R 



A quantity which describes how the estimates J2k=i ^jkSk of Yj couple to the 
measurements Yi is given by the na x na matrix 



a(EE|M^ ,,,, 

dYi dY ^ ' 

also called regularization matrix [17]. The trace of R can be interpreted as the 
effective number of measurements used in the unfolding procedure. The num- 
ber of effectively rejected measurements is then equal to Nj-ej — na—Trace{R). 



The maximum value of Trace{R) , which is equal to the rank nr of G, is reached 
with the direct solutions (4) and (7), corresponding to the cases nr — na < nb 
and na > nb — nr respectively: For the direct solution 5"° the number of 
rejected measurements Nrej is equal to na — Trace{R) = na — nr — na- 
na = 0, which means that no information is discarded. The measurements Y 
are completely reproduced by the unfolding : M-S^ = Y. For the least squares 
solution the number of effectively rejected measurements is Nrej ~ 

na — Trace{R) — na — nr — na — nb > 0, which means that some 
information is discarded. The measurements Y are not exactly reproduced by 
the unfolding : M-S^^^ ^ Y . The fact that the system is overconstrained has 
a similar effect as regularization. In both cases, with increasing regularization 
strength Trace{R) is reduced and Nrej is increased. 



6 Selecting the unfolding result 



For a given unfolding method the result S depends on the regularization 
strength, which is given by the weight w (or by the number of iterations i 
respectively). In the literature various criteria for choosing the "best" weight 
are proposed [8,10,13,16-18]. Unfortunately, none of them provides a choice 
which is optimal for all cases. Reasons for this are: The optimum regulariza- 
tion strength in general depends on the shape of the unknown distribution S. 
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It also depends on the binnings in Etme and Eggt and on the prior distribution 
(if appUcable) . 

The effect of the regularization is illustrated in Fig. 1, where different quan- 
tities are plotted as a function of the iteration number i. In this example an 
experimental energy distribution of 7-rays from the Crab Nebula [24], i.e. the 
number of excess events in bins of the estimated energy E^st, was unfolded. 
The unfolding was performed for 30 different i in the range 10~^ to 10^°, using 
the method of Bertero (eq.(20)). More results from the analysis of these data 
are given in Section 10. With decreasing i, i.c with increasing regularization 
strength, one observes an increase of Xo (^1- 3) and a decrease in the quantities 
Trace{T) /Trace{K), Trace{R), Reg{S)Tikhonov (eq. H) and Reg{S) schmeiung 
(eq. 13). Very similar behavior is found for the other unfolding methods, dis- 
cussed in Section 3. 

Obviously, an acceptable unfolding result should satisfy the following condi- 
tions: 

• The x^-probability, calculated from the value of Xq and the number of de- 
grees of freedom in the unfolding procedure, should be acceptable, say > 1%. 
Otherwise the unfolding result is incompatible with the measured distribu- 
tion Y . 

• The noise term Trace{T) of the unfolded distribution S should be compa- 
rable to the noise term Trace{K) of the measurements. The main aim of 
regularization is a suppression of the large noise term of S, which one often 
obtains if no regularization is applied. A large noise term Trace{T), as well 
as large correlation terms of T, indicate a too fine binning in Etrue, leading 
to small eigenvalues of G or if (see Section 7.1). 

• Trace{R) should not be much lower than its maximum possible value, which 
is equal to the rank nr of the matrix G. Otherwise the solution is too strongly 
dominated and biased by the regularization. 

For determining the "best" regularization strength a compromise has to be 
found between the above requirements. It has turned out that the criterion 
Trace{T) = Trace{K) in general leads to solutions which satisfy the above 
conditions reasonably well, provided the problem is not strongly overcon- 
strained. In the latter case, where the unfolding result is better constrained, 
a solution with Trace{T) < Trace{K) is more apropriate. In MAGIC 
the standard criterion for determining the optimal regularization strength is 
Trace(T) = Trace{K). The full circles in Fig. 1 indicate this choice. However, 
any other regularization strength can be chosen by hand, if this is suggested 
by the behaviour of the quantities Xq, Trace{T), Trace{R), Reg{S)Tikhonov 
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Fig. 1. Useful quantities for determining the optimum regularization strength, plot- 
ted as a function of the iteration number i. 
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or Reg{S) Schmelling- 



Unfolding with regularization is a procedure which allows freedom in the choice 
of the regularization method and in the choice of the regularization strength. 
The above criteria for an acceptable solution strongly restrict this freedom. 
Nevertheless a certain degree of arbitrariness remains as to which unfolding 
result should be considered representative and final. In MAGIC a selected 
unfolding result is considered representative if all other unfolding methods 
yield results, which are also acceptable and statistically consistent with the 
selected result. In addition, it is required that also the Forward Unfolding 
(Section 4), using a reasonable parametrization of ^{Etrueji gives a consistent 
result. 

An uncertainty due to the unfolding is determined from the spread of the Sj, 
obtained from the different unfolding methods. 



7 Further comments on the unfolding 

In the actual application of the unfolding procedure to real data some technical 
details have to be considered, and they are discussed in this Section. 



7.1 Optimal binnings 

The binning of the experimental distribution Y is often dictated by the avail- 
able statistics and by the experimental errors. The binning should not be 
chosen too fine in order to assure significant measurements in all bins. In 
the case of an lACT experiment, a sufficiently large sample is required to 
determine the number of signal (excess) events with sufficient accuracy. The 
binning should not be chosen too wide either, because the binning in Y limits 
the reconstruction of the fine structure of the unfolded distribution S. 

Another criterion for the binnings is the behavior of Gram's matrix G. A too 
fine binning for S leads to strong correlations between neighboring columns 
of the migration matrix M, implying small eigenvalues of G, which lead to a 
large noise component of S. Given a certain choice of the binning in E^gf and 
thus of na, the bin size in Efrue or nb should be set such, that the system of 
linear equations (2) is not underconstrained. This usually leads to wider bins 
in Etrue than in Eest- In MAGIC a typical value of Alogio{Etrue) / ^iogio{Eest) 
is 1.4 (see Section 10). 
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It should also be noted that the unfolding procedure doesn't require equidis- 
tant bins, neither in Eggt nor in Et^ue- 



7.2 Completeness of the migration matrix 

If nal < i < na2 specifies the range of bins of the measured distribution 
Y which are to be considered in the unfolding procedure, also the range in 
Etrue to be considered in the unfolding has to be chosen properly: one has to 
make sure that all bins j of Ef^ue ^-re present, for which the column j of the 
migration matrix contributes to the selected bins of Y, i.e. for which at 
least one of the elements Mij {nal <i< na2) is different from zero. 

There is an exception to this rule, if for some reason certain bins j of E^rue 
are not expected to contribute to the selected bins in Eggt- This is for example 
the case if one of the factors A^^^ in (22) is so small that Sk can be neglected. 



7.3 Additional cuts 

As explained in Section 2.1, the distribution Y and the migration matrix 
M have to be produced under identical conditions (selections, cuts). If an 
additional cut is imposed when generating Y, also M has to be recalculated 
before doing the unfolding. If this additional cut is a cut in Egst one may 
proceed in the following way: 

• Renormalize the columns j of M to the selected range in E^st (see eq.(l)). 

• Perform the unfolding of Y in the usual way. 

• Apply a correction A-'^^^^^ to the unfolded distribution ^j, where A-'^^^^^ is 
the renormalization factor for column j of M. 



7.^ Starting values for the minimization 

In the cases where the unfolding procedure involves numerical minimizations, 
like the ones discussed in section 3.1, the minimization may not converge. 
This problem can be often solved by choosing different starting values. Another 
reason for non-convergence is discussed in Section 7.5. In MAGIC the standard 
choice of the starting distribution for S*, and also of the prior distribution 
e in Schmelling's method, is a distribution which is close to the measured 
distribution Y . 
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7. 5 Components of S which cannot be determined in the unfolding 



According to cq. (2) those Sj for which the column Mij (?' = 1, ...na) is a 
null vector have no influence on Y and can therefore not be determined in 
the unfolding. These components should not be varied in the minimization 
because they would unnecessarily complicate the minimization process and 
may lead to non-convergence. 



7.6 Dependence on the assumptions made in the Monte Carlo simulation 



For the unfolding the migration matrix Mij is the crucial quantity. Obviously, 
if it doesn't describe the real migration of events correctly, the unfolding result 
will be wrong. This means that at flxed j, i.e. at flxed Efrue, the MC simulation 
has to describe the migration in E correctly. This will be the case if at fixed 
Etrue the shower simulation is realistic and if the detector response is simulated 
correctly. 

On the other hand, the distribution of Etme in the MC need not agree with the 
real distribution of E : due to the normalization of M (eq. (1)) the bin-to-bin- 
variation of the number of MC events in Efrue has no influence on Mij at all. 
This is one of the great advantages of unfolding methods like those presented 
in Section 3 as compared to methods based on correction factors (see Section 
11). 

However, there is a residual dependence of Mij due to the flnite binning in 
Etrue '■ Depending on the shape of the Etme distribution within an Etme bin 
in the MC simulation, the calculated Mij may be more representative for the 
lower, middle or upper part of the Etme bin. If the Etme distribution in the 
real data is different from that in the MC simulation the calculated M^j may 
not be exactly the right ones. 

This residual dependence of M on the shape of the Etme distribution in the 

MC can be nearly completely removed by an iteration procedure in which the 
M for the next iteration step is determined from a MC sample, in which the 
Etme distribution has been corrected using the unfolding result of the last 
iteration step (see Section 8). 
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8 Determining the effective collection cirea A and the migration 
matrix M for a finite bin in Et^ue ^ind © 



In an lACT experiment it is important that the effective collection area A, 
which enters in the ffux calculation, is computed talcing into account a realistic 
shape of the differential flux ^{Etme) and the actual distribution of effective 
observation times dT{Q)/dQ in the zenith angle 0. This is also important for 
the migration matrix M, which enters in the unfolding, because of the residual 
dependence on the flux spectrum, as discussed in Section 7.6. Recalculating A 
and M with the proper Efrue and spectra is the more important the bigger 
the (AEtrue, A0) interval and the stronger the variations of M and A within 
this bin are. Often, for statistics reasons, large bin sizes in Eest and 0, and 
thus also in Etme, have to be chosen. 

In the following it is assumed that the effective collection area A and the 
migration matrix M are known functions of Efrue and 0. This can be achieved 
by determining them from a sample of MC 7-ray events in very flne bins of 
Etrue and 0. The aim is to calculate an average A of the effective collection 
area and an average M of the migration matrix, which are representative for 
a flnite bin {AEtrue, A©) in Etrue and 0. 



8.1 The effective collection eirea 



The number of observed events in a {AEtrue, A©) bin is given by 

S ^ [ f AiEtrue, 0) ■ HEtrue) ' " dE^me " dQ (25) 

JAG J AEtrue CLKy 

Here ^(Etrue) is the differential 7-ray ffux to be measured, is the 

distribution of observation times in the experimental data and A(Etrue, 0) is 
the known dependence of the effective collection area on Etme and 0. 

With the deflnitions 

AT = / ^^^^ ■ dO (= total observation time) (26) 

^AO 00 



$ = -TT; / HEtrue) ■ dEtrue (27) 

J AEt rue 



Ia0 i AEtrue ^i-^true, ©) " ^(Etrue) " J^q ^ " dEf^ue " dQ 



AT • $ ■ AE, 



(28) 



true 
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equation (25) can be rewritten as 



$ = = (29) 

AT-A-AEtrue 

This is the usual formula for converting numbers of events S into differential 
fluxes Because of the definitions (27) and (28) the differential flux $ in an 
Etrue bin, as determined in the unfolding, is the average differential flux in 
this bin. Therefore, when quoting or plotting a result for <l> the bin edges in 
Etrue should also be given or shown. 



8.2 The migration matrix 

The number of reconstructed MC events in a bin i of Eest can be written as 

Ni ^ f f MiiEtrue, e) -AiEtrue, 6) ■ ^Etme) ' ^^^^ " dEtme " dQ 

J AO JAEtrue 

(30) 

with ^Mk{Etrue,Q) = 1 for all E^true and © 

k 



Mi{Etrue,&) is the element of the normalized migration matrix for the i-th 
bin in Egst , at an energy Efme- The index j of Mij is replaced by the variable 
Etrue- The dependence of Mi{Etrue, ©) on Etrue and © is assumed to be known 
from MC simulations. 

The average migration matrix Mj for the selected AEtrue bin, to be used in 
the unfolding, is then obtained by 

_ N, 
Mi = 



EkNk 



Ia0 I AEtrue ^ii^true, ©) " M^true, ©) " ^{Etrue) ' " dEtme ' dQ 



IaO JAEtrue ^i^true, ©) " ^{Etrue) ' " dEtrue ' dO 



(31) 



The averages and A are calculated according to (31) and (28) respectively, 
using an approximation ^i{Etrue) of the function ^ (Etrue)- The measured dis- 
tribution Y is unfolded using Mj, yielding the unfolded distribution 5". A new 

approximation $2(-E'true) is then determined from S according to eq. (29). The 
procedure is iterated until ^(Etrue) has converged. In practice parametric func- 
tions are used as approximations of ^{Etrue), and in general the convergence 
is found to be very fast. 
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9 Combining data before applying the unfolding procedure 



One often has the situation that there exist several measured distributions 
Y" of the same quantity. If the Y" were obtained under different conditions 
also the migration matrices M" and the effective collection areas A" will be 
different for the different measurements. 

In an lACT experiment the different conditions may be 

• Different modes of observation (ON/OFF mode, wobble mode, observation 
in the presence of moon light, ...). 

• Different ranges of the zenith angle. 

• Different detector conditions. 

• etc. 

In order to determine a final unfolded distribution S one may proceed in 
different ways : 

• Individual unfolding : 

Unfold each Y'^ using M" to obtain S'^, and combine the S" to obtain the 
final solution S. 

• Global unfolding : 

Combine the Y'^ and to obtain a global Y and M, do an unfolding of 
Y using M, which will give the final solution S. 

There is one important argument in favour of the second option : Each of the 
Y" , or some of them, may have large statistical errors, making the unfolding 
of the individual Y"^ unstable. One common unfolding of the global Y using 
the global M will be more robust, in general. 

In the case of an I ACT experiment, the following relations hold for each mea- 
surement u : 



j denotes the j-th bin in Etrue, is its width, is the average flux and 
Aj is the effective collection area in this bin, and T is the effective observation 
time. Since $ and are the same for all v the relations for the combined 
data read : 



Y 




T" ■ A"^ ■ ■ AEj 



(32) 



Y = M-S 



(33) 
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Inserting Y ^J^i^Y" and T ^ ^^T" in (33) and using (32) one obtains 



A, - ^^^^ (34) 

^ _ E.M-<-(A-n _ E.M-.-(A-n 

The relations (34) and (35) give the prescription how to combine the individual 
to obtain the global M, and how to combine the individual Aj to obtain 
the global Aj. The measured distribution Y is unfolded using M, and the 
unfolded distribution S is converted into a flux $ according to (33). As can be 
seen from (34) and (35) the Aj are weighted averages of the Aj with weights 
= T'^ /T, and the Mij are weighted averages of the M^j with weights 
Wi, — A^jT^ I {J^i^ A^T^^). The weights can be interpreted as fractions because 
they add up to 1. They are correlated and their covariance matrix has to be 
taken into account when calculating the errors of Aj and M^j. 

The equations (34) and (35) also show that Aj and M^j can be obtained with- 
out knowning which spectra ^{Etrue) and dT{Q)/dQ were used to compute 
the A", and M'!-. 



10 Application to experimental data 



In this Section an experimental energy distribution of 7-rays from the Crab 
Nebula is unfolded, which was obtained in an analysis of data taken with the 
MAGIC telescope [24]. The migration matrix M, as determined from a sample 
of 7-MC events, is plotted in Fig. 2a) as a function of Eest and Etme- The size 
of the boxes is proportional to the value of Mij, where i and j are the bin 
numbers in E^gt and Efrue respectively. The experimental distribution Y of 
the number of 7-excess events as a function of E^st is displayed in Fig. 2b). 
Both distributions are after all cuts and selections, except a cut in Eest- 

The vertical and horizonthal lines in Figs. 2a) and c) indicate the ranges 
in Eest or Efrue, which were selected for the unfolding. The range in Eest is 
given by those Eest bins for which a significant number of excess events could 
be determined. The range in Etme comprises all those Etme bins which are 
expected to contribute to the selected range in Eest- Fig- 2a) suggests that these 
are all 15 Etme bins. However, according to the plot of the effective collection 
area Aejj in Fig. 2c) the contribution from Etrue < 60 GeV is expected to be 
neghgible. The same holds for Etrue > 9 TeV, due to the strongly decreasing 
7-ray flux with increasing energy. This is confirmed by a Forward Unfolding of 

dN ( E \" 

Y in which a differential 7-ray fiux of the form — ; — = fo • I ^ — I , 

' ^ dA-dt-dE ■'^ V300 GeVy ' 
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Fig. 2. Input data for the Unfolding and for the calculation of the differential 7-ray 
flux. 



with a = a + 6 ■ lo£ 



E 



was assumed. This leads to a number of 



V300 GeV, 

Eest and Etme bins of na = 14 and nb = 10 respectively, which are used in the 
unfolding procedure. The size of the logio{Etrue) bins was deliberately chosen 
wider than the logio{Eest) bins by a factor of 1.4, in order to better constrain 
the unfolding. The rank of Gram's matrix G is equal to nr = 10, as can be seen 
from Fig. 2d), which shows the size of the eigenvalues A; of G as a function 
of /. Two to three of these eigenvalues are much smaller than the maximum 
eigenvalue, and they are the reason for the large values of Trace{T) /Trace{K) 
at large iteration number (low regularization strength) in Fig. lb). 

The optimum regularization strength and thus the final solution S was deter- 
mined using the criterion TraceiT) /Trace{K) = 1. The estimates J2j ^ij " Sj 
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(a) Unfolded distribution S before (grey (b) Result for the differential 7-ray flux 
symbols) and after the correction for the ^{Etme) multiplied with E"^. 
cut in Eest (black symbols). 
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(c) Comparison of the results for E'^ ■ 
^(Etrue) from different unfolding meth- 
ods. 



Fig. 3. Results from the Unfolding. 



(open circles) of Yi are compared with Yi (histogram) in Fig. 2b). These data 
enter in the calculation of Xo (^Q- 3). The number of degrees of freedom in the 
unfolding procedure is na — nb + 1 = 5, because the number of measurements 
is na, the number of unknowns is nb, and the relation J2j Sj = J2i Yi (eq. 6) is 
used as additional constraint. As can be seen from Figs. Id) and e) the values 
of Reg {S) Tikhonov and Reg{S)schmeUing are much lower at the selected regular- 
ization strength than without regularization. This means that the solution 5* 
is smoothed by the regularization. The black symbols in Fig. 3a) represent the 
final solution S. The solution before correctiong for the cut in i^gst is drawn 
with grey symbols. 
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The final differential 7-ray flux as computed from S according to eq.(22), is 
drawn in Fig. 3b). The solid line represents the result of a fit of the expression 

300 Gev ) ' ^^^^ a = « + ^' ■ logio ( 300 Cev ) ' ^^^^ points. The 

number of degrees of freedom for this fit is n6 — 3 = 7, because the number 
of data points is nb = 10 and the number of free parameters is 3 (/o, a and 
b) . The is 8 for 7 degrees of freedom. Setting 6 = in the fit yields a of 
24 for 8 degrees of freedom. This fit is clearly disfavoured as compared to the 
fit in which the slope a is energy dependent. In these fits the full correlation 
matrix T oi S has been taken into account. 

The result of the latter fit for ^{Efrue) was used to recalculate the averages 
Mi (31) and A (28) for the individual Efrue bins. The unfolding was repeated 
using the recalculated Mi and A, yielding new results for S, $ and the fit 
parameters. After 1 iteration this procedure converged. 

Very similar results were obtained with the other unfolding methods. The 
spread of the Sj, obtained with the different unfolding methods, can be seen 
in Fig. 3c). This spread can be regarded as an estimate of the a systematic 
error due to the unfolding. 



11 Correction Factors 



A widely used method of correcting experimental distributions is the applica- 
tion of correction factors: Using Monte Carlo data, both the 'true' distribution 
S^'-^ {k = 1, ...nc) and the 'reconstructed' distribution Yj^'-^ {k = 1, ...nc) 
of some quantity are produced under certain conditions (selections, cuts). Cor- 
rection factors are determined according to 

Ck = S^IY^^ {k = 1, ...nc) (36) 

An experimental distribution Y^, obtained under the same conditions as , 
is then corrected by 



Sk ^ Yk- Ck {k = 1, ...nc) (37) 

to obtain the corrected distribution Sk- 

The following properties of this procedure can be stated [17] : 
Ck is undefined if Y,^'^ = 0. 



If Sj^^ — also Ck — and Sk — 0, which means that is ignored. 
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• Cfc depends on the shape of the MC distribution S^'-^; the corrected distri- 
bution Sk is always biased towards S^'-^. 

• If Yfc = also Sk is zero. 

• The standard linear error propagation often yields too small errors of Sk- 

The correction factors are only right if Sj^'^ is identical to the true Sk distri- 
bution. If this is not the case one may iterate S^'" , setting S^^ equal to the 
last corrected experimental distribution Sk- However, this often leads to insta- 
bilities. The reason for the instabilities appears to be similar to that causing 
a large noise component of the direct solution (4) . 

In contrast to the unfolding methods presented in Section 3, there is very little 
freedom in choosing the binnings for 5* and Y- By definition, the range of values 
and the binnings for the true and reconstructed quantity are identical. 

Advantages of the method of correction factors are that it is simple and stable. 
The drawbacks have been listed above, the severest one being the strong de- 
pendence of the correction factors on the assumptions made in the MC about 
S- 



12 Summciry 



In this paper the procedures to unfold experimental energy distributions of 7- 
rays, as applied in the MAGIC experiment, are described. It is explained, how 
the uncertainties, which are inherent in any unfolding process, can be handled 
successfully. Possible problems in the unfolding are discussed and suggestions 
are given which can help to avoid them. Various techniques are presented, 
which allow to reconstruct the energy spectrum in a rather unbiased way. All 
algorithms are impleneted in the MAGIC software, which is based on the (7"^+ 
language and ROOT [25]. Their application to real data has shown to provide 
robust and reliable results. The methods and procedures are applied in most 
of the MAGIC analyses. 
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